PenmanMonteith Subroutine

SUBROUTINE PenmanMonteith &
!
(airTemp, netRad, rh, wind, fc, z, elevation, zws, zrhum, lai, srmin, pt, pe, pet)

IMPLICIT NONE

!Arguments with intent(in):
REAL (KIND = float), INTENT(in)  :: airTemp !!air temperature [°C]
REAL (KIND = float), INTENT(in)  :: netRad !! net radiawion [W/m2]
REAL (KIND = float), INTENT(in)  :: rh !! air relative humidity [0-100]
REAL (KIND = float), INTENT(in)  :: wind !! wind speed [m/s]
REAL (KIND = float), INTENT(in)  :: fc !! fractional coverage by vegetation [0-1]
REAL (KIND = float), INTENT(in)  :: z !! vegetation height [m]
REAL (KIND = float), INTENT(in)  :: elevation !! terrain elevation [m a.s.l.]
REAL (KIND = float), INTENT(in)  :: zws !! wind speed measurement heigth [m]
REAL (KIND = float), INTENT(in)  :: zrhum !! relative humidity measurement heigth [m]
REAL (KIND = float), INTENT(in)  :: lai !! leaf area index [m2/m2]
REAL (KIND = float), INTENT(in)  :: srmin !! minimum stomatal resistance

!Arguments with intent(out):
REAL (KIND = float), INTENT(out) :: pt		!! potential transpiration (from vegetation) [m/s]
REAL (KIND = float), INTENT(out) :: pe		!! potential evaporation (from water or saturated soil) [m/s]
REAL (KIND = float), INTENT(out) :: pet		!! potential evapotranspiration [m/s]

!local declarations:
REAL (KIND = float) :: groundHeat !!ground heat flux [MJ m-2 s-1]
REAL (KIND = float) :: es !!saturation vapor pressure (Pa)
REAL (KIND = float) :: ea !!actual vapor pressure (Pa)
REAL (KIND = float) :: des !!the slope of the relationship between saturation vapour pressure and temperature (kPa°C-1)
REAL (KIND = float) :: airPress !!air pressure [KPa]
REAL (KIND = float) :: gamma !!psychrometric constant [kPa °C-1]
REAL (KIND = float) :: ws !!wind speed (m/s)
REAL (KIND = float) :: z1  !!bare soil heigth (m)
REAL (KIND = float) :: zm1  !! wind speed measurement heigth + z1 (m)
REAL (KIND = float) :: zh1  !! relative humidity measurement heigth + z1 (m)
REAL (KIND = float) :: zd1   !! bare soil zero plane displacement height (m)
REAL (KIND = float) :: zrw1  !! Roughness length governing momentum transfer above zd1 (m)
REAL (KIND = float) :: zrh1   !! Roughness length governing transfer of heat and vapor above zd1 (m)
REAL (KIND = float) :: rabs !! aerodynamic resistance of bare soil (s/m)
REAL (KIND = float) :: zm   !! wind speed measurement heigth + z (m)
REAL (KIND = float) :: zd    !! vegetated soil zero plane displacement height (m)
REAL (KIND = float) :: zrw   !! Roughness length governing momentum transfer above zd (m)
REAL (KIND = float) :: zh  !! relative humidity measurement heigth + z (m)
REAL (KIND = float) :: zrh   !! Roughness length governing transfer of heat and vapor above zd (m)
REAL (KIND = float) :: ra !! aerodynamic resistance of vegetated soil (s/m)
REAL (KIND = float) :: cp !!  humid air specific heat (MJkg^-1°C^-1)
REAL (KIND = float) :: Tkv  !! virtual air temperature (K)
REAL (KIND = float) :: rhoAir !!air density (kgm^-3)
REAL (KIND = float) :: netRadMJ !! net radiation in MJm^-2s^-1
REAL (KIND = float) :: rc !! total vegetation  resistance (s/m)
REAL (KIND = float), PARAMETER :: k = 0.41 !!von Karman’s constant
REAL (KIND = float), PARAMETER :: epsilon = 0.622   ! ratio molecular weight of water vapour/dry air
REAL (KIND = float), PARAMETER :: lambda = 2.453 !2.2647 !!latent heat of vaporization (MJ / kg)
REAL (KIND = float), PARAMETER :: R1 = 0.287	!!gas specific constant	(kJkg^-1K^-1)
!-----------------------------------------end of declarations------------------

! compute saturation vapor pressure
es = 0.6108 * EXP ( ( 17.27 * airTemp ) / ( airTemp + 237.3 ) )

!actual vapor pressure
ea = rh / 100. * es

!compute the slope of saturation vapour pressure curve (Allen. et al. 1998)
des = ( 4098. * es ) / ( ( airTemp + 237.3 )**2. )

!compute atmospheric pressure (ASCE-EWRI, 2005)
airPress = 101.3 * ( ( 293. - 0.0065 * elevation ) / 293. )** 5.26

!compute psychrometric constant
gamma = 0.665 * 0.001 * airPress

!check wind speed: set minimum to 0.01 m/s
IF ( wind <= 0 ) THEN
	ws = 0.01
ELSE
	ws = wind
END IF

!aerodynamic variables of bare soil
z1 = 0.1 
zm1 = zws + z1
zd1 = 2. / 3. * z1
zrw1 = 0.123 * z1
zh1 = zrhum + z1
zrh1 = 0.1 * zrw1

!compute aerodynamic resistance of bare soil
rabs = LOG ( (zm1 - zd1) / zrw1 ) * LOG ( (zh1 - zd1) / zrh1 ) / ( k**2 * ws ) 

!aerodynamic variables of vegetated soil
zm = zws + z
zd = 2. / 3. * z
zrw = 0.123 * z
zh = zrhum + z
zrh = 0.1 * zrw

!compute aerodynamic resistance of bare soil
IF (z == 0 ) THEN  !vegetation not present
	ra = rabs
ELSE
	ra = LOG ( ( zm - zd ) / zrw ) * LOG ( ( zh - zd ) / zrh ) / ( k**2 * ws)
END IF

! humid air specific heat (MJkg^-1°C^-1)
cp = gamma * epsilon * lambda / airPress

!virtual air temperature (K)
Tkv = 1.01 * ( airTemp + 273. )

!air density (kgm^-3)
rhoAir = airPress / ( Tkv * R1 )


!check and convert net radiation
IF ( netRad < 0 ) THEN
	netRadMJ = 0.
ELSE
	netRadMJ = netRad * 0.000001
END IF 

! compute ground heat flux as a fraction of net radiation Anderson et al. (2007)
! groundHeat is assumed to affect evaporation from bare soil
!groundHeat = 0.31 * ( 1. - fc ) * netRadMJ
groundHeat = 0.1 * netRadMJ 


!compute potential evaporation (from saturated bare soil or water)
pe = ( des * ( netRadMJ - groundHeat ) + rhoAir * cp * ( es - ea ) / rabs  )  / ( lambda *  ( des + gamma )  ) * millimeter
    

!compute potential transpiration from vegetation
IF (fc > 0.) THEN  
    IF ( lai == 0.) THEN
        rc = srmin / 1.
    ELSE
        rc = srmin / lai
    END IF
    
    pt = ( des * netRadMJ + rhoAir * cp * ( es - ea ) / ra  )  / ( lambda *  ( des + gamma * (1. + rc / ra) )  ) * millimeter

ELSE
    
    pt = 0.

END IF

!compute potential as a weighted average of bare soil and vegetated fluxes
pet = ( 1. - fc) * pe + fc * pt

!IF (isnan (pet)) then
!    write(*,*) pe, pt, des, ra, rc, lai, srmin
!    pause
!end if

RETURN
END SUBROUTINE PenmanMonteith